function check (type)
% Check variables
%
global nlp
global jpi hpi

if strcmp(type, 'jacobian') == 1
    checkJac();
elseif strcmp(type, 'hessian') == 1
    checkHessian();
end
%% functions
    function checkJac ()
        xJ = nlp.guess; rj = randperm(nlp.nvar, 10);
        xJ(rj) = 1.1 * xJ(rj);
        % --- Jacobian by "analytic"               --- %
        J = feval(nlp.func.jacobian, xJ);
        % --- Jacobian by "finite-differentiation" --- %
        c = feval(nlp.func.constraints, xJ);
        Jf = zeros(nlp.ncon, nlp.nvar);
        turb = 10*sqrt(eps);
        for i = 1 : nlp.nvar
            xi = xJ; xi(i) = xi(i) + turb;
            ci = feval(nlp.func.constraints, xi);
            Jf(:,i) = (ci-c)/turb;
        end
        Jf = sparse(Jf);
        fprintf('Max Jacobian error of J(analytic) - J(finite-diff):\n');
        fprintf('    %f\n', max(full(max(J-Jf))));
        % --- Jacobian by "auto-differentiation"   --- %
        if isempty(jpi) == 1
            jpi = getjpi(nlp.func.constraints, nlp.nvar);
        end
        [~, Ja] = evalj(nlp.func.constraints, xJ, [], nlp.ncon, jpi);
        fprintf('Max Jacobian error of J(analytic) - J(auto-diff):\n');
        fprintf('    %f\n', max(full(max(J-Ja))));
    end
    function checkHessian ()
        lmbd = ones(nlp.ncon, 1); rl = randperm(nlp.ncon, 10);
        lmbd(rl) = 1.1 * lmbd(rl);
        xh = nlp.guess; rx = randperm(nlp.nvar, 10);
        xh(rx) = 1.1 * xh(rx);
        % ---       Hessian by "analytic"         --- %
        H = feval(nlp.func.hessian, xh, [], lmbd);
        % --- Hessian by "auto-differentiation"   --- %
        if isempty(hpi) == 1
            hpi = gethpi('lagAD', nlp.nvar);
        end
        [~, ~, Ha] = evalh('lagAD', xh, {'', lmbd}, hpi);
        Ha = tril(Ha);
        fprintf('Max Hessian error of H(analytic) - H(auto-diff):\n');
        fprintf('    %f\n', max(full(max(H-Ha))));
    end
end